Thermodynamic properties and structural stability of thorium dioxide 
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Using density functional theory (DFT) calculations, we have systematically investigated the ther- 
modynamic properties and structural stabilities of thorium dioxide (TI1O2). Based on the calculated 
phonon dispersion curves, we calculate the thermal expansion coefficient, bulk modulus, and heat 
capacities at different temperatures for TI1O2 under the quasi-harmonic approximation. All the 
results are in good agreement with corresponding experiments proving the validity of our methods. 
Our theoretical studies can help people more clearly understand the thermodynamic behaviors of 
TI1O2 at different temperatures. In addition, we have also studied possible defect formations and 
diffusion behaviors of helium in TI1O2, to discuss its structural stability, ft is found that in intrinsic 
TI1O2 without any Fermi energy shifts, the interstitial Th* + defect other than oxygen or thorium va- 
cancies, interstitial oxygen, and any kinds of Frenkel pairs, is most probable to form with an energy 
release of 1.74 eV. However, after upshifting the Fermi energy, the formation of the other defects 
also becomes possible. For helium diffusion, we find that only through the thorium vacancy can it 
happen with the small energy barrier of 0.52 eV. Otherwise, helium atoms can hardly incorporate 
or diffuse in TI1O2. Our results indicate that people should prevent upshifts of the Fermi energy of 
TI1O2 to avoid the formation of thorium vacancies and so as to prevent helium caused damages. 

PACS numbers: 63.20.dk, 65.40.-b, 61.72.-y, 66.30.J-. 
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I. INTRODUCTION 



As the world's demands for energy keep growing, cor- 
responding researches on developing new energy sources 
or enhancing the energy-consuming efficiency are attract- 
ing more and more attentions. Until now, the fissile nu- 
clear reactor is still a very important energy source, in 
which uranium dioxide (UOV) has been the main fuel 
component for many years |l|-_5|. However, during the 
burning cycle of UO2, considerable amounts of redioac- 
tive elements emerge in the reaction waste @ . And such 
radioactive waste results in very troublesome long-term 
storage requirements. Based on these facts, many ef- 
forts have been done to look for possible substitutions 
of UOV Thorium based materials, which are naturally 
abundant, are identified as good candidates for replacing 
UO2 in fissile nuclear reactors because they are able to 
produce fewer transuranic (TRU) compared to uranium- 
and plutonium-based fuels. And although 233 Th is not 
fissile, it can absorb a slow neutron and then form fissile 
233 U just undergoing multiple beta-decays [7j]. Moreover, 
because of the good solubility between thorium dioxide 
(Th02) and other transuranic dioxides, using the mix 
oxides (MOX) of thorium and plutonium in nuclear reac- 
tors can also help reducing the large plutonium stockpile 
while maintaining acceptable safety and control charac- 
teristics of the reactor system [8[. In comparison with 
previous uranium-based fuels, the thorium-based fuels 
also have many additional physical advantages, such as 
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higher melting points, higher corrosion resistivity, lower 
thermal expansion coefficients, and higher thermal con- 
ductivity [9J. In recent years, the thorium-based fuels 
have already been tested in different reactors @, . 

As in the front position of actinide element, thorium 
and its compounds have been studied ever since 1950. At 
the earliest stage, a series of experimental measurements 
have been carried out on the thermal properties of Th02 
such as the thermal expansion coefficient, heat capacity, 
and thermal conductivity fUBW-j-j . In 1997, Bakker et al. 
(23| presented a conclusion and comparison for the mea- 
sured values. On the other side, Chadwick and Graham 
[23 |, Allen and Tucker [Hj|, and Veal et al. [H] investi- 
gated the valence-band structures of thorium and its ox- 
ides by means of X-ray photoemission spectroscopy. And 
the pressure-induced phase transition of TI1O2 has been 
studied by Jayaraman [27} , Dancausse [28[ , and Idiri [29[ 
et al. experimently. Recently in 2006, researchers from 
the international atomic energy agency (IAEA) built a 
thermal-physical database of materials for light water re- 
actors and heavy water reactors [30j . And the earlier 
experimental measurements are re-addressed. 

Despite the vast experimental measurements, it is to 
our surprise that no one has ever theoretically investi- 
gated the thermodynamic properties and defect behav- 
iors for TI1O2, which are critically important for its usage 
in thermonuclear reactors. Only recently, several theo- 
retical studies have been carried out on the mechanical 
and electronic properties [3ll - [33| . phase transition be- 
haviors [13, and elastic and optical properties [36[ 
for TI1O2. Our previous study has already proven that 
the thorium 5/ states is no longer localized after elec- 
tronic hybridizations, and density functional theory cal- 
culations are enough to produce correct descriptions for 
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the ground-state properties of TI1O2 [35(. So in our 
present paper, we decide to systematically investigate 
the thermodynamic properties and structural stabilities 
of ThC>2, by using density functional theory calculations. 
The thermodynamic stability will be discussed based on 
the calculated thermal parameters, while the structural 
stability will be discussed by calculating the formation 
energy of different kinds of defects, and diffusion energy 
barriers of helium in ThC>2. The rest of the paper is 
organized as follows. The computation methods are in- 
troduced in Section II. The discussions about the thermo- 
dynamic properties and structural stabilities of TI1O2 are 
presented in Section III. Finally, we give our conclusions 
in Section IV. 



II. CALCULATION METHODS 

The density functional theory (DFT) calculations are 
carried out usin g th e Vienna ah initio simulations pack- 
age (vasp) [al |38[ with the projector-augmented-wave 
(PAW) potential methods [39J . The cutoff energy for the 
plane-wave basis set is set to 500 eV. The exchange and 
correlation effects are described by generalized gradient 
approximation (GGA) in the Perdew-Burke-Ernzerhof 
(PBE) form ■ A 2 x 2 x 2 supercell is employed to study 
defect formation and helium diffusions inside TI1O2. For 
calculations of the unit cell (12 total atoms) and 2x2x2 
supercell (96 total atoms), the integration over the Bril- 
louin Zone is done on 13x13x13 and 5x5x5 fc-ponit 
meshes generated using the Monkhorst-Pack [4l[ method, 
which are both proven to be sufficient for energy conver- 
gence of less than l.OxlO -4 eV per atom. During the 
supercell calculations, the shape and size of the super- 
cell are fixed while all the ions are free to relax until the 
forces on them are less than 0.01 eV/A. 

For a semiconductor, the Hclmholtz free energy F at 
volume V and temperature T can be expressed as 



defined as, 



F{V,T) = E{V) + F vib {V,T), 



(1) 



where E(V) stands for the ground-state electronic en- 
ergy, F v n,(V,T) is the phonon free energy at a given 
unit cell volume V. Within quasi-hamonic approxima- 
tion (QHA), F V ib(V,T) can be evaluated by 



F vib {V,T) 



3,1 



2 sinh 



V 2fc B T 



(2) 



where u)j (q, V) is the phonon frequency of the jth phonon 
mode with wave vector q at fixed V, and k B is the Boltz- 
mann constant. The total specific heat of the crystal is 
the sum of all phonon modes over the Brillouin zone (BZ), 



C v (T) = J2cvA<hT). 



(3) 



c^(q,T) = k B ^2 



3,q 



2k B T 



sinh 2 [/^ (q, V)/2k B T] 
(4) 

The mode Gruneisen parameter 7j(q) describing the 
phonon frequency shift with respect to the volume can 
be calculated by 

#n Wj (q,VO] 



7j(q) 



(5) 



d[Ln V] 

The acoustic Gruneisen parameter j(T) defined as the 
weighted average of the mode Gruneisen parameters for 
all acoustic phonon branches is calculated to be 



a v (T)B(T)V m (T) 



(6) 



Cu.j(q, T) is the mode contribution to the specific heat 



C V (T) 

where a v (T) is the thermal expansion coefficient and 
equals to y (§jr)p) V m (T) is the volume per mole ma- 
terial, and B(T) and C V (T) are the bulk modulus and 
specific heat respectively. 

The formation energy of a point defect X with charges 
q can be calculated by introducing the chemical potential 
concept as, 

E for {X") = E tot ±n x n-E Th02 +q{E v + E f + AV), (7) 

where E to t is the total energies of the supercell with 
defect X, n x represents the number of X defects, [i is 
the chemical potential of X with a positive sign for va- 
cancy and a negative sign for interstitial defect, i?Tho 2 
is the energy of the TI1O2 supercell without defects, E v 
and Ef are the valence-band maximum (VBM) and the 
Fermi level of Th02 respectively. With these denota- 
tions, E for (V o ), % or (09), £ /or (V° h ), and £/or(Th°) 
represent for the formation energies of a neutral oxygen 
vacancy, a neutral interstitial oxygen, a neutral thorium 
vacancy, and a neutral interstitial thorium respectively. 
The shift of the VBM in a defect supercell AV takes the 
change of the valence-band maximum caused by the de- 
fect into account. Its value can be obtained by a macro- 
scopic average technique [HI, |43[ through calculating the 
difference between the average electrostatic potential in 
a bulklike environment of the defect supercell and the av- 
erage electrostatic potential in the defect-free supercell. 
The formation energy of a Frenkel-pair can be calculated 
by summing up the formation energies of a vacancy and 
a corresponding interstitial add- in, i.e., 

E for (FP X ) = E for (V X ) + £/cr(Xi). (8) 

For the Schottky defect, the formation energy can be 
calculated by 

Efor(S) — S/or(Vxh) + 2Ef or(Vo) — — ^-E T h0 2 , 

(9) 

where N is the number of atoms in the considered super- 
cell. In this expression, the defect consists of a thorium 
vacancy and two oxygen vacancies, which are again sup- 
posed to be non-interacting. 
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III. RESULTS AND DISCUSSIONS 
A. Structure and Elastic constants of TI1O2 

Our previous studies have shown that DFT calcu- 
lations with GGA exchange-correlation Junctionals are 
good enough for obtaining the ground-state properties of 
ThC>2 [47], and adding additional U or J modifications 
for localization effects might lead to incorrect results. We 
believe that it is because the 5/ electronic states of Th be- 
come delocalized after electronic hybridizations in TI1O2. 
Here based on our previous studies, we further investigate 
the thermodynamic properties of TI1O2. 

From the mechanical point of view, the theoretical 
equilibrium volume Vb, bulk modulus Bq, and the pres- 
sure derivative of bulk modulus B' can be obtained by 
fitting the third-order Brich-Murnaghan equation of state 
[44| . In this way, our calculated lattice parameter a for 
ThC>2 is 5.619 A, which is in accordance with the exper- 
imental data of 5.60 A [H, Sj|. The bulk modulus B Q 
and its pressure derivative B' are calculated to be 190 
GPa and 4.3, also in agreement with correspon ding ex- 
perimental values of 195-198 GPa and 4.6-5.4 (HSg, re- 
spectively. In order to evaluate the Poisson's ratio v : we 
calculated the three independent elastic constants Cn, 
C%2 and C44 of Th02- The calculation methods are the 
same as in our previous studies [H, EH, S3- The ob- 
tained elastic constants for TI1O2 are Cn=351.2 GPa, 
Ci2=106.9 GPa, and 644=74.1 GPa, which are in accor- 
dance with experimentally measured values of Cn=367 
GPa, Ci2=106 GPa, and C* 44 = 79 GPa 0. Furthermore, 
the Poisson's ratio v is calculated to be 0.293, in excel- 
lent accordance with the experimental data of 0.285 [Hj]. 
All the good agreement between our calculated values 
and corresponding experimental measurements indicates 
that our calculation methods are effective and reliable. In 
comparison with other actinide dioxides, we can see that 
Th02 has a slightly smaller bulk modulus than UO2 and 
Pu0 2 (471. 
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FIG. 1: (Color online) Phonon dispersion curves and phonon 
density of states (DOS) for Th0 2 . 
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FIG. 2: (Color online) Temperature dependence of the linear 
thermal expansion for TI1O2. The inset is the volume thermal 
expansion coefficient as a function of temperature. 



B. Phonon dispersions, Thermal expansion, and 
Heat capacity of TI1O2 

The thermodynamic properties of a material are con- 
nected to its phonon dispersion curves. In order to get 
accurate splittings between longitudinal (LO) and trans- 
verse optical (TO) phonon branches, the Born effective 
charge is firstly calculated. Because of the high symme- 
try of Th02, the Born transverse effective charge ten- 
sor Z* is the same along the [100], [010], and [001] di- 
rections, and the effective charge can be averaged by 
Z* = \TrZ*. For Th0 2 , we obtain that Z* h =+5.41, 
and Zq=— 2.71. The static dielectric constant e = \Tre 
is 4.83. The dielectric constant might be smaller than 
experimental measurements because of the underestima- 
tion of the electronic energy band gap due to the draw- 
back of the exchange-correlation approximation (GGA). 



The calculated phonon dispersion curves along the high- 
symmetry k-point lines using the above Born effective 
charges are shown in Fig. 1 . As clearly shown, there 
is an obvious splitting between longitudinal optical (LO) 
and TO branches due to the Born effective charge. More- 
over, we find that there is no evident gap between acous- 
tic and optical branches of phonon for Th02 , with an ob- 
servable overlap between the longitudinal acoustic (LA) 
and transverse optical (TO) branches around the X point. 
Detailed vibrational modes analysis tells us that the vi- 
brations of Th and O atoms dominate the low- (0-6 THz) 
and high-frequency (6-18 THz) modes, respectively. 

From the obtained phonon dispersion curves and Eqs. 
(1) and (2), we then calculate the energy curves for TI1O2, 
and find the lowest-energy lattice constants at different 
temperatures. The lattice expansion curve is thus ob- 
tained and shown in Fig. 2. The experimental data by 
Touloukian et al. [49[ and Kim et al. [30| are also shown 



4 



o 

60 

~3 

-a 
o 

m 



198 



192 



186 



180 



174 



168 



1 1 1 1 1 1 


I I 


1.00 






— 0.96 

- e 

_ =0.92 






~~ 0.88 
- 






300 600 900 1200 1500 
T(K) 

I , , I , , I , , I , , 



300 600 900 1200 1500 
Temperature (K) 



100 



i4 80 



o 
E 



60 



40 



20 



_ ■ ■ ■ ■ 1 ■ 

- 


i i i | i i i i | i i i i | i i i i ^ 
i 


- 


^oSSSSSS— —— — : 

Cv -This work - 




Cp -This work . 




■ Bakker et al. - 




• Tim et al. 


i 


" 



300 600 900 
Temperature (K) 



1200 1500 



FIG. 3: (Color online) Temperature dependence of bulk mod- 
ulus B for TI1O2. The inset is the ratio of B /Bo. 



FIG. 4: (Color online) Heat capacities at constant volume 
(C„) and constant pressure (C p ) of Th02- 



in Fig. 2 for comparisons with our calculational results. 
One can see that in the temperature range from 300 to 
600 K, our result is in excellent agreement with the ex- 
perimental data. In a higher temperature range from 
600 to 1500 K, our obtained theoretical values are in the 
middle of the two different experimental measurements. 
And at 1500 K, the relative difference between our result 
and the two experimental values are 0.15% and 0.05%, 
respectively. The small relative differences indicate that 
the QHA method can give reasonable lattice parameters 
for Th02 up to 1500 K. The thermal expansion coeffi- 
cient a v (T) is calculated and shown in the inset of Fig. 
2. The experimental values for the thermal expansion co- 
efficient in the temperature range from 298 to 1500 K can 
be calculated from the corresponding thermal expansion 
curves, which are 3.318 xl0 -5 i4T -1 for Touloukian's data 
and 3.630 xlO~ 5 K^ 1 for Kim's data, respectively. Our 
theoretical result of 3.509 xlO^AT -1 is consistent with 
the experimental values. 

The bulk modulus B is also analyzed as a function of 
temperature according to the formula B — Vo(^y?)v , 
and the result is displayed in Fig. 3. We can see that the 
value of bulk modulus decreases with increasing temper- 
ature, and at 1500 K, the ratio of bulk modulus (B/Bq) 
is 0.898, as depicted in the inset of Fig. 3. 

Under QHA, the considered vibration modes are har- 
monic but volume-dependent. The calculated heat ca- 
pacity at constant volume using Eqs. © and Q is shown 
in Fig. 4, together with the heat capacity at constant 
pressure C p , which is calculated according to the rela- 
tionship, 



C p - C v = a 2 v (T)B(T)V(T)T. 



(10) 



The available experimental data from Bakker et al. [23| 
and Kim et al. [30j are also shown in Fig. 4 for compar- 
isons. It is evident that our theoretical result is in excel- 
lent agreement with the measured values for the whole 
experimentally considered temperature range. As tem- 
perature increases, the value of C p increases continuously, 



while the value of C v approaches to a constant of 3R (R 
is the Rydberg constant). At 1500 K, the value of C p be- 
comes 81 J-mol _1 K _1 . In general, QHA method is valid 
when the temperature is much lower than the material's 
melting point (around 3600 K for TI1.O2 [23j) when the 
anharmonic effect is small. 



Griineisen parameters and thermal conductivity 
of Th0 2 



The lattice thermal conductivity kl for a material can 
be calculated differently, depending on the specific mech- 
anisms for phonon scattering. At relative high temper- 
atures, the dominant mechanism for phonon scattering 
is the Umklapp process, in which the acoustic phonon 
branches interact with each other and transport heat. 
With this mechanism, the lattice thermal conductivity 
of a crystal-like material can be expressed as [5bTj52| , 



K L =A 



Me 3 (T)5(T)n 2 / 3 
7 2 (T) x T ' 



(11) 



where A is a physical constant with the value of 
3.1xl0~ 6 , M is the average mass per atom in the crys- 
tal, 0(T) is the Debye temperature of TI1O2, n is the 
number of atoms in the primitive unit cell, "f(T) is the 
acoustic Griineisen parameter, and S(T) is the cube root 
of the average volume per atom, i.e., the averaged radius 
per atom. The kl and 6 in Eq. ([TTjl are in units of W 
m — 1 K _1 , and A respectively. With reasonable expres- 
sions of the Debye temperature and acoustic Griineisen 
parameter to describe the harmonic phonon branches and 
the anharmonic interactions between different phonon 
branches, Eq. ([lip can provide accurate predictions for 
a material's thermal conductivity [5l[ . 

Firstly the Debye temperature O can be determined 
from the elastic constants within the Debye theory, in 
which the vibrations of the solid are considered as elastic 
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FIG. 5: (Color online) The thermal conductivity of TI1O2. 
Exp erimental results from Springer et al. Weibacher 
McElroy et al. Ql], ARF El, Murabayashi et al. [3, 
Pillai and Raj Ql, Bradshaw [2(j, Moore et al. [2l[, Murti 
and Mathews [22||, and Bakker et al. [23J are displayed for 
comparison. 



waves, and the Debye temperature of the solid is related 
to an averaged sound velocity [Hj]. Within isotropic ap- 
proximation, the Debye temperature can be expressed 
as 15311 . 



0(D = ^ V ^ { T)n]^f{ V )^l^, (12) 

where M is the molecular mass per formula unit, B(T) 
is the bulk modulus, v is the material's Poisson's ratio, 
and f{v) is given by 0, El] 



/0) = 3 1/3 



2 1 



3 1-2^ 



3/2 



1 1 



3 1 - v 



3/2' 



-1/3 



(13) 

Our calculations find that the Debye temperature mono- 
tonically decreases, while the acoustic Griineisen parame- 
ter monotonically increases with increasing temperature. 

Based on Eq. (fTTj) and the obtained Debye temper- 
ature, we calculate the thermal conductivity for Th02- 
Figure 5 shows our thermal conductivity result, in com- 
parison with the previous experimental measurements in 
the temperature range from 300 to 1500 K. We can see 
that in the considered temperature range, our calculated 
values are in agreement with the experimental results, 
which proves the validity of our methods and model. Es- 
pecially, in the high temperature range from 600 to 1200 
K, our results accords very well with the experimental 
measurements by Murabayashi, et al. [l8j and by Pil- 
lai and Raj [l9j]. At the low temperature range around 
the Debye temperature of 402 K, our results are slightly 
different from some experimental results. This difference 
comes from our presumption that the dominant mecha- 
nism for phonon scattering is the Umklapp process. The 



FIG. 6: (Color online) The mode Griineisen parameters along 
high-symmetry directions in the reciprocal lattice space of 
TI1O2. 



accordance between our calculations and corresponding 
experiments in the temperature range from 600 to 1200 
K proves that within this temperature area, the contri- 
butions from other phonon-scattering mechanisms are so 
small that can be neglected. 

The mode Griineisen parameter describing the phonon 
frequency shift with respect to the volume can be used 
to discuss the anharmonic effects. By expanding or com- 
pressing the equilibrium volume by 1%, we calculate the 
mode Griineisen parameter 7j(q) for all nine phonon 
branches according to Eq. ([5]). The corresponding re- 
sults are shown in Fig. 6. It can be seen that all the 
mode Griineisen parameter values are positive, indicat- 
ing that all phonon frequencies increase with decreasing 
volume. Besides, the acoustic phonon mode Griineisen 
parameters are relatively larger reflecting that changes 
in volume have more influences on the collective vibra- 
tion modes of TI1O2. We can also see from Fig. 6 that the 
LA and TO phonon branches have larger mode Griineisen 
parameters, indicating that the anharmonic interactions 
between the LA and TO branches should be more inten- 
sive with respect to the volume change. 



D. Defect formation in ThC>2 

In this subsection and the next, we will investigate the 
structural stability of Th02 by systematically calculating 
the formation energy of different kinds of defects and 
investigating the diffusion behaviors of helium. A 2 x 2 x 2 
supercell with 96 Th and O atoms is employed in these 
two subsections to model defect formation and helium 
diffusion in Th02 . Different charge states are considered 
for all the defects. To calculate the formation energy, the 
positions of all ions are fully relaxed before we calculate 
the electronic free energies of the system with different 
defects. 

The formation energy for oxygen (thorium) vacan- 
cies, interstitial oxygen (thorium) ions, oxygen (thorium) 
Frenkel pairs, and the Schottky defect of a TI1O2 unit in 
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TABLE I: Formation energies of different defects at different charge states in TI1O2. The defects include oxygen vacancy (Vo), 
interstitial oxygen ion (O;), thorium vacancy (VtiO, interstitial thorium ion (Thi), oxygen (FPo) and thorium Frenkel-pairs 
(FPxh), and Schottky defect (S) of a TI1O2 unit. The formation energies are in units of eV. 
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FIG. 7: (Color online) Formation energies of oxygen related defects, i.e., oxygen vacancies or interstitial oxygen ions (a) and 
thorium related defects, i.e., thorium vacancies or interstitial thorium ions (b) as functions of the Fermi energy. 



the 2x2x2 supercell are obtained and listed in Table fl] 
One can see that the charged Thf* + add-in defect is the 
only one having a negative formation energy. It means 
that once Th* + ions are available, interstitial Th* + de- 
fect can form in TI1O2 with an energy release of 1.74 eV. 



Comparatively for oxygen defects, the O 2- vacancy is 
the most possible one because of its smallest formation 
energy among oxygen defects. These two results reflects 
the ionic character of the Th-0 bonds that each tho- 
rium atom almost loses 4 electrons to two oxygen atoms 
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FIG. 8: (Color online) Electronic density of states (DOS) for TI1O2 with oxygen vacancies (a), and with interstitial oxygen 
defects (b) in different charge states. The Fermi energies are set to zero. 



in Th02- From Table HI we can also see that all tho- 
rium vacancies and thorium Frenkel-pairs have huge for- 
mation energies and thus can hardly form in intrinsic 
ThC>2 • Comparatively, the formation of oxygen Frenkel- 
pair, containing an 2+ vacancy and a interstitial O 2- 
ion has a relatively smaller formation energy. 

Depending on the environment's influence on the Fermi 
energy of the Th02 material, the formation possibility of 
charged defects can change due to the changes on hard- 
ness of electrons transfer. Using Eqs. (7)- (9) under differ- 
ent Fermi energies, we can calculate the formation energy 
as a function of the Fermi energy for different kinds of 
defects. The corresponding results for oxygen (thorium) 
vacancies and interstitial oxygen (thorium) ions in differ- 
ent charge states are respectively shown in Figs. 7(a) and 
7(b). The experimentally measured band-gap of 6.0 eV, 
instead of the 4.1 eV value obtained by first-principles 
calculations is chosen as the reference band-gap in dis- 
cussions on defect formation. 

As can be seen from Fig. 7(a), when comparing the en- 
ergies associated with forming V Q , Vq, and Vq + oxygen 
vacancies, a transition can be observed with the Fermi 
level increasing from the valence band to the conduction 
band. The +2 charged oxygen vacancy is favored near 
the valence band, indicating that oxygen vacancies have 
a tendency to donate electrons or behave as a n-type de- 
fect. When the Fermi level increases to around 2.6 eV, 
the Vq becomes energetically favorable. With further in- 



creasing the Fermi level, the neutral charge state is most 
probable for oxygen vacancies, and the tendency of oxy- 
gen vacancy to donate electrons diminishes. We can also 
see from Fig. 7(a) that the interstitial 2 ~ defect can 
become very possible to form when the Fermi energy is 
shifted to be near the conduction band of TI1O2. 

From Fig. 7(b) we see that when the Fermi energy 
is near the valence band, all charged states of thorium 
vacancy are hard to form because of the huge formation 
energies, and the interstitial Th^ + ions can easily form in 
TI1O2 with a negative formation energy. With upshift- 
ing the Fermi energy, the interstitial thorium defect in 
neutral state and the V^Jj vacancy can both become the 
most possible thorium kinds of defects. When consid- 
ering only the vacancy defects, V Q + is the most stable 
defect near the valence band, while near the conduction 
band, V^ is the most favorable one. Comparatively for 
interstitial states, the Th^ + and 2 ~ defects are the most 
probable ones when the Fermi energy is near the valence 
and conduction band respectively. 

In addition, the electronic density of states (DOS) for 
both defect-free and defective TI1O2 are calculated and 
shown in Figs. 8 and 9, to further analyze the influences 
of the considered defects on the electronic structures of 
Th02- As clearly shown in both Figs. 8 and 9, the in- 
troduction of vacancies or interstitial ions do not change 
the DOS distribution of TI1O2 very much. The biggest 
character in both Figs. 8 and 9 is that a new defect en- 



8 



160 - 

80 - 

- 
160 - 

80 - 

- 
160 - 

80 - 

1 



w 160 - 
O 

a 

80 - 

- 

160 - 

80 - 

- 

160 - 

80 - 



Defect-free 



-20 -16 -12 -8 -4 4 
(a) Energy (eV) 



160 
80 


160 
80 



Defect-free 



Th° 



oi^ 



160 - 

? 80 - 
c 

3 

■e 



Th. + 







g 160 
Q 

80 



TV 





160 - 
80 - 


160 - 
80 - 



Th. 







-20 -16 -12 -8 -4 4 
(b) Energy (eV) 



FIG. 9: (Color online) Electronic density of states (DOS) for TI1O2 with thorium vacancies (a), and with interstitial thorium 
defects (b) in different charge states. The Fermi energies are set to zero. 



ergy level emerges in the band gap for defective TI1O2. 
For the electronic structure of TI1O2 with oxygen vacan- 
cies, we see from Fig. 8(a) that the Fermi energy shifts 
from the valence band maximum (VBM) in defect-free 
ThC>2, to above the defect energy level for the Vq defect, 
and from above the defect energy level for the Vq defect 
back to the VBM for the Vq + defect. For the system 
with interstitial oxygen ions, we can see from Fig. 8(b) 
that the Fermi energy shifts respectively from the VBM 
for the 0° defect to the defect energy level for the 0~ 
defect, and from the defect energy level for the de- 
fect to above the defect energy level for the O^ - defect. 
Similarly for the TI1O2 supercell with interstitial thorium 
ions, the Fermi energy shifts from the VBM to the defect 
energy level for the Th+, Th?+, and Th?+ defects, and 
then shifts to above the defect energy level for the Th? + 
defect, as shown in Fig. 9(b). From Fig. 9(a), one can 



see that the defect energy levels for thorium vacancies in 
TI1O2 are very close to the VBM. 



E. Diffusion of helium in TI1O2 

In order to determine the structure influence of helium 
impurity on Th02, we systematically calculate the incor- 
poration energy and diffusion energy barriers for helium 
in both the defect-free and defective TI1O2. The incor- 
poration energy is defined to be the energy required to 
incorporate one helium atom at a pre-existing vacancy or 
at an interstitial site. In this way, E inc can be expressed 
as follows 



E, 



tot 



^Th0 2 — ^He, 



(14) 
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FIG. 10: (Color online) Diffusion pathways for a helium atom from one octahedral interstitial site to another octahedral 
interstitial site directly in intrinsic TI1O2 (a), through an oxygen vacancy (b) and through a thorium vacancy in defective ThC>2 
(c). 
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FIG. 11: (Color online) The energy profiles for a helium atom to diffuse from one octahedral interstitial site to another 
octahedral interstitial site directly in intrinsic TI1O2 (a), through an oxygen vacancy (b) and through a thorium vacancy in 
defective ThQ 2 (c). 



where E to t is the energy of the TI1O2 supercell with an 
incorporated helium, i?Th0 2 is the energy of the Th02 su- 
percell with or without defects, and En e the energy of an 
isolated helium atom. Three different kinds of incorpo- 
rations are systematically considered here, i.e., a helium 
atom at the oxygen (O) and thorium vacancy (Th) sites 
in defective Th02, and a helium atom at the octahedral 
interstitial site (OIS). The calculated incorporation ener- 
gies are 2.48, 0.21, and 0.77 eV respectively. We can see 
that the most possible incorporation site for helium is the 
thorium vacancy site in defective Th02- In comparison, 
the oxygen vacancy site is much less possible for helium 
incorporation. 

The structural influence of helium on Th02 is deter- 
mined not only by the incorporation energy, but also 
by the minimum energy path for helium diffusion. The 
climbing image nudged elastic band (CINEB) method 
[56| is then employed to find the minimum-energy dif- 
fusion pathways. Based on the obtained incorporation 
results for helium in TI1O2, we here investigate the dif- 



fusions of helium from OIS to OIS in both defect-free 
and defective TI1O2. Figures 10 (a)- (c) depict our consid- 
ered diffusion ways. And the energy profiles along these 
paths are shown in Figs. ll(a)-(c) respectively. We can 
see from Fig. 11(a) that in intrinsic TI1O2, the incorpo- 
rated helium atom needs to overcome a 3.34 eV energy 
barrier in order to diffuse from one OIS to another. This 
huge diffusion energy barrier means that at normal tem- 
peratures, helium diffusion in intrinsic Th02 is almost 
impossible. In the corresponding saddle point along this 
diffusion path, the helium atom is at the middle of the 
two OISs, and the nearest oxygen atoms along the [001] 
direction are repelled from each other by about 0.80 A. 

From the calculated energy profile shown in Fig. 11(b), 
we see that the energy barrier for an incorporated helium 
atom to diffuse from one OIS to another by passing an 
oxygen vacancy is lowered down to be 1.58 eV. This en- 
ergy barrier is however still too large for the diffusion to 
happen at room temperatures. The only possible way for 
helium to diffuse in Th02 is found to be passing through 
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a thorium vacancy, as depicted in Fig. 10(c). We can 
see from Fig. 11(c) that the helium atom only needs 
to overcome a 0.10 eV energy barrier to diffuse from 
an OIS to the thorium vacancy, and a 0.52 eV energy 
barrier to diffuse from the thorium vacancy to another 
OIS. The small energy barriers for helium to diffuse in 
thorium vacancy-included TI1O2 indicate that thorium 
vacancies might lead to helium aggregation causing fail- 
ure of the TI1O2 material. Fortunately from our above 
defect formation studies, the formation of thorium va- 
cancies is almost forbidden when the Fermi energy is not 
upshifted, because of their too large formation energy in 
ThC>2. Therefore, to keep TI1O2 away from structural 
damages from helium incorporation, any kinds of factors 
possibly leading to upshifts of the Fermi energy should 
be avoided. 



IV. SUMMARY 

In summary, we have performed a systematic first- 
principles study to investigate the thermodynamic prop- 
erties and structural stabilities of TI1O2. Based on the 
calculated phonon dispersion curves for TI1O2, we sys- 
tematically analyze its thermodynamic properties and 
obtain the values of its thermal expansion coefficient, 
bulk modulus, and heat capacities at different temper- 
atures, which are in good agreement with correspond- 
ing experimental measurements. The agreement between 
our calculations and experiments also proves the validity 
of our methods and model, and the effectiveness of the 
quasi-harmonic approximation. According to the Umk- 
lapp interaction mechanism between different phonon 



branches, we systematically obtain the mode Griineisen 
parameters, and further calculate the thermal conductiv- 
ities of TI1O2. Within the temperature range from De- 
bye temperature to about 1500 K, our calculated thermal 
conductivity accords very well with experimental results. 

In addition to studying the thermodynamic properties 
of ThC>2, we also investigate its structural stability by 
calculating the formation energy of different defects, and 
the diffusion behaviors of helium, during which different 
charge states of the defects are considered. The forma- 
tion energy results indicate that without any shifts of the 
Fermi energy, the interstitial Th 4+ defect is very prob- 
able to appear in TI1O2 with an energy release of 1.74 
eV. With changing the Fermi energy to different values, 
the formation possibilities of different defects varies. For 
helium incorporation, it is found that the helium atom 
tends to occupy a thorium vacancy in defective TI1O2 or 
occupy the octahedral interstitial site in intrinsic Th02- 
It is further revealed that incorporated helium atoms 
can only diffuse freely in the thorium-vacancy contained 
ThC>2, with small energy barriers of 0.10 and 0.52 eV. 
Our studies point out that to avoid helium damage, the 
electronic Fermi energy of TI1O2 should not be upshifted 
because it can makes the formation of thorium vacancies 
less possible. 
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